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QQ ' Abstract 

o ' 

In Part I of this paper we discussed new methods for the numerical continuation of point-to- 
cycle connecting orbits in 3-diniensional autonomous ODE's using projection boundary condi- 
^ ' tions. In this second part we extend the method to the numerical continuation of cycle-to-cycle 

j^ ■ connecting orbits. In our approach, the projection boundary conditions near the cycles are 

formulated using eigenfunctions of the associated adjoint variational equations, avoiding costly 
and numerically unstable computations of the monodromy matrices. The equations for the 
eigenfunctions are included in the defining boundary-value problem, allowing a straightforward 
implementation in AUTO, in which only the standard features of the software are employed. 
Homotopy methods to find the connecting orbits are discussed in general and illustrated with 
an example from population dynamics. Complete AUTO demos, which can be easily adapted to 
any autonomous 3-dimensional ODE system, are freely available. 

Keywords: boundary value problems, projection boundary conditions, cycle-to-cycle connec- 
tions, global bifurcations. 



1 Introduction 

In a diversity of scientific fields bifurcation the- 
ory is used for the analysis of systems of ordi- 
nary differential equations (ODE's) under param- 
eter variation. Many interesting phenomena in 
ODE systems are linked to global bifurcations. 
Examples of such are overharvesting in ecologi- 
cal models with bistability properties (Bazykin, 
1998; Antonovsky et al., 1990; Van Voorn et al., 
2007), and the occurrence and disappearance of 
chaotic behaviour in such models. For example, 
it has been shown (see Kuznetsov et al., 2001 and 
Boer et al., 1999, 2001) that chaotic behaviour 
of the classical food chain models is associated 
with global bifurcations of point-to-point, point- 
to-cycle, and cycle-to-cycle connecting orbits. 

In Part I of this paper (Doedel et al., 2007) we 
discussed heterodinic connections between equi- 
libria and cycles. Here we look at connections 
that link a cycle to itself (a homodinic cycle-to- 
cycle connection, for which the cycle is necessarily 
saddle), or to another cycle (a heterodinic cycle- 
to-cycle connection). Orbits homoclinic to the 
same hyperbolic cycle are classical objects of the 
Dynamical Systems Theory. It is known thanks 
to Poincare (1879), Birkhoff (1935), Smale (1963), 
Neimark (1967), and L.P. Shilnikov (1967) that a 
transversal intersection of the stable and unsta- 
ble invariant manifolds of the cycle along such an 
orbit implies the existence of infinite number of 
saddle cycles nearby. Disappearance of the in- 
tersection via collision of two homoclinic orbits 
[homodinic tangency) is an important global bi- 
furcation for which the famous Henon map turns 
to be a model Poincare mapping (Gavrilov and 
Shilnikov, 1972; Palis and Takens, 1993, see also 
Kuznetsov, 2004). 

Numerical methods for homoclinic orbits to 
equilibria have been devised by Doedel and Kern- 
evez (1986, but see Doedel et al., 1997), who ap- 
proximated homoclinic orbits by periodic orbits 
of large but fixed period. Beyn (1990) developed 
a direct numerical method for the computation 
of such connecting orbits and their associated pa- 



rameter values, based on integral conditions and 
a truncated boundary value problem (BVP) with 
projection boundary conditions. 

The continuation of homoclinic connections in 
AUTO (Doedel et al., 1997) improved with the de- 
velopment of HomCont by Champneys and Kuz- 
netsov (1994) and Champneys et al. (1996). How- 
ever, it is only suited for the continuation of bifur- 
cations of homoclinic point-to-point connections 
and some heterodinic point-to-point connections. 
A modification of this software was introduced by 
Demmel et al. (2000), that uses the continuation 
of invariant subspaces (CIS-algorithm) for the lo- 
cation and continuation of homoclinic point-to- 
point connections. 

Died and Rebaza (2004) have also made sig- 
nificant progress recently, by developing methods 
to continue point-to-cycle and cycle-to-cycle con- 
necting orbits based on another work by Beyn 
(1994). Their method employs a multiple shoot- 
ing technique and requires the numerical solving 
for the monodromy matrices associated with the 
periodic cycles involved in the connection. 

Our previous paper (Doedel et al., 2007) dealt 
with a method for the detection and continuation 
of point-to-cycle connections. Here this method 
is adapted for the continuation of homoclinic and 
heterodinic cycle-to-cycle connections. The me- 
thod is set up such that the homoclinic case is 
essentially a heterodinic case where the same pe- 
riodic orbit (but not the periodic solution) is dou- 
bled. In Section [2] we give a short overview of 
a BVP formulation to solve a heterodinic cycle- 
to-cycle problem. In Section [3] it is shown how 
boundary conditions are implemented. In Sec- 
tion H] we discuss starting strategies to obtain ap- 
proximate connecting orbits using homotopy. In 
Section [5] the BVP is made suitable for numerical 
implement at ion . 

Results are presented of the continuation of a 
homoclinic cycle-to- cycle connection in the stan- 
dard three-level food chain model in Section El 
Boer et al. (1999) previously numerically obtained 
the two-parameter continuation curve of this con- 
necting orbit using a shooting method, combined 



with the Poincare map technique. In the previ- 
ous part of this paper (Doedel et ah, 2007) we 
reproduced the results for the structurally stable 
heteroclinic point-to-cycle connection of the same 
food chain model using the homotopy method. In 
this paper we discuss how the homoclinic cycle-to- 
cycle connection can be detected, and continued 
in parameter space using the homotopy method. 
Also, it is set up such that it can be used as well 
for a heteroclinic cycle-to-cycle connection. 

2 Truncated BVPs with pro- 
jection BCs 

Before presenting the BVP that describes a cycle- 
to-cycle connection, let us first set up some nota- 
tion. Consider a general system of ODEs 

f = /K«), (1) 

where / : R" x M^ — > M" is sufficiently smooth, 
given that state variables u G M", and control 
parameters a G MP. Thus, the dimension of the 
state space is n and the dimension of the param- 
eter space is p. The (local) flow generated by ([1]) 
is denoted by </?*. Whenever possible, we will not 
indicate explicitly the dependence of various ob- 
jects on parameters. 

We assume that both 0~ and O^ are saddle 
limit cycles of ([1]). A solution u{t) of ([1]) for fixed 
a defines a connecting orbit from 0~ to O''" if 

lim dist(M(t),0^) =0 . (2) 

t— »itoo 

(Figure [T] depicts such a connecting orbit in the 
3D-space.) Since u(t + r) satisfies ([T]) and ([2]) for 
any phase shift r, an additional phase condition 



Table 1: List of notation used in the paper. 



sym. 


meaning 


x± 


Periodic solution 


v^ 


Eigenf unction 


w^ 


Scaled adjoint eigenfunction 


u 


Connection 


a 


Bifurcation parameters 


0+ 


Limit cycle where connection ends 


0- 


Limit cycle where connection starts 


wi 


Stable manifold of the cycle 0+ 


w 


Unstable manifold of the cycle 0~ 


f^t 


Unstable multiplier of the cycle 0~^ 


f^s 


Stable multiplier of the cycle 0~ 


f^u 


Unstable multiplier of the cycle 0~ 


f^^ 


Adjoint multiplier l//i^ 


/i 


Adjoint multiplier l//i7 


A± 


ln(/.±) 


T± 


Period of the cycle 0^ 


T 


Connection time 



is prescribed that the points u{t^) and u{t+) be- 
long to the linear subspaces, which are tangent to 
the unstable and stable invariant manifolds of 0~ 
and O"*", respectively. 

Now, denote by x^ (t) a periodic solution (with 
minimal period T^) corresponding to O^ and in- 
troduce the monodromy matrix 



M^ 



Dx^ 



T±, 



[X 



= (0) 



I.e. 



■ip[u, a] = 



(3) 



should be imposed to ensure uniqueness of the 
connecting solution. This condition will be spec- 
ified later. 

For numerical approximations, the asymptotic 
conditions ([2D are substituted by projection bound- 
ary conditions at the end-points of a large trun- 
cation interval [t-,T-^], following Beyn (1994). It 



the linearization matrix of the T^'^-shift along 
orbits of ([T]) at point Xq = x^(0) G O^. Its eigen- 
values are called Floquet multipliers, of which one 
(trivial) equals 1. Let mf = nj" -t- 1 be the di- 
mension of the stable invariant manifold W^ of 
the cycle 0~^, where n^ is the number of its mul- 
tipliers satisfying 

|/i| < 1. 

Along the same line, m~ = n~ -|- 1 is the dimen- 
sion of the unstable invariant manifold W^ of the 
cycle 0~, and n~ is the number of its multipliers 



W'^/ 



o 





(a) 



(b) 



Figure 1: Cycle-to-cycle connecting orbits in 
orbit, 0+ = 0-. 



p3. 



heterolinic orbit, O^ i^ O ; (b) honioclinic 



satisfying 



|/i| > 1. 



To have an isolated branch of cycle-to-cycle 
connecting orbits of ([1]) it is necessary that 



V 



n — m' 



m. 



(4) 



(see Beyn, 1994). 

The projection boundary conditions in this 
case become 



L^{u{t^)-x^{Q)) =0 



(5) 



where L~ is a (n — m^) x n matrix whose rows 
form a basis in the orthogonal complement to the 
linear subspace that is tangent to W^ at x"(0). 
Similarly, L"*" is a (n — rn^) x n matrix, such that 
its rows form a basis in the orthogonal comple- 
ment to the linear subspace that is tangent to 
W'^ at x+(0). 

The above construction also applies in the case 
when 0+ and 0~ coincide, i.e. we deal with a 
homoclinic orbit to a saddle limit cycle 0~^ = 0~. 
Note that, in general, the base points x^(0) G O^ 
remain different (and so do the periodic solutions 



X 



± 



m 



It can be proved that, generically, the trun- 
cated BVP composed of ([1]), a truncation of ([3]), 
and ([5]), has a unique solution branch {u{t, a), a), 



provided that ([T]) has a connecting solution branch 
satisfying ([3]) and (jl]). 

The truncation to the finite interval [r_ , r^] 
causes an error. If u is a generic connecting so- 
lution to ([1]) at parameter a, then the following 
estimate holds in both cases: 



(u 



.^+1' 



a) 



(u, a) 



< (^g-2min(^_|r-.|,^+|T+|) 



(6) 

where || ■ || is an appropriate norm in the space 
C*^([r_,r+],]R") X Rp, ti|[r_,r+] is the restriction of 
u to the truncation interval, and n± are deter- 
mined by the eigenvalues of the monodromy ma- 
trices. For exact formulations, proofs, and refer- 
ences to earlier contributions, see Pampel (2001) 
and Dieci and Rebaza (2004, including Erratum). 

3 New defining systems in M^ 

In this section we show how to implement the 
boundary conditions ([5]). We consider the case 
n = 3 where 0~ and O"*" are saddle cycles and 
therefore always m~ = ?n~ = 2 and mf = m^ = 
2. Substitution in (jlj) gives the number of free 
parameters for the continuation p = 1. 

Note that the complete BVP will consist of 15 
equations (2 saddle cycles, 2 eigendata for these 
cycles, and the connecting orbit) and 19 boundary 
conditions. 




Figure 2: Ingredients of a BVP to approximate a heteroclinic connecting orbit. The homoclinic cycle- 
to-cycle connection is also approached as the heteroclinic case, where two saddle cycles coincide. 



3.1 The cycle and eigenf unctions 

To compute the saddle limit cycles 0~ and O^ 
involved in the heteroclinic connection (see Fig- 
ure [2]) we need a BVP. The standard periodic 
BVP can be used 



X- 



-fix' 



,a 



a;±(0)-a;±(T±) 








(7) 



A unique solution of this BVP is determined by 
using an appropriate phase condition, which is 
actually a boundary condition for the truncated 
connecting solution, and which will be introduced 
below. 

To set up the projection boundary condition 
for the truncated connecting solution u near O^, 
we also need a vector, say w~^{0), that is orthog- 
onal at x~^{0) to the stable manifold W^ of the 
saddle limit cycle O"*", as well another vector, say 
w~{0), that is orthogonal at a;~(0) to the unsta- 
ble manifold PV"" of the saddle limit cycle 0~ (see 
Figure [2]). Each vector w^{0) can be obtained 
from an eigenf unction w'^it) of the adjoint varia- 
tional problem associated with (jTj), corresponding 
to eigenvalue /x^. These eigenvalues satisfy 



/i^ 



K 



/^ 



where //^ and /i^ are the multipliers of the nion- 
odromy matrix M^ with 



The corresponding BVP is 

'^ + Pu 



< 1 . 



^Tfx±, 



w±(T±)-/i±w±(0) 
(u7^(0),u;^(0)) 



1 









(8) 



where x"^ is the solution of ([7]). In our implemen- 
tation the above BVP is replaced by an equivalent 
BVP 

+ f'^{x^,a)w^ + \^w^ 

{w^{0),w^{0))-l 



w"^ 









(9) 



where s 



± 



sign fi"^ and 



A^ 



f^'s 



ln|/i±| . 

(See Appendix of Part I, Doedel et al., 2007). 

In Q), the boundary conditions become peri- 
odic or anti-periodic, depending on the sign of the 
multiplier /x^, while the logarithm of its absolute 
value appears in the variational equation. This 
ensures high numerical robustness. 

Given w"^ satisfies ([9]), the projection bound- 
ary conditions ([S]) become 

(w;±(0),M(r±)-a;±(0)) = 0. (10) 



3.2 The connection 

We use the following BVP for the connecting so- 
lution: 



ii — f{u, a) 
(/(x±(0),a),n(r±)-x±(0)) 



0, 
0. 



:iii 



For each cycle, a phase condition is needed to se- 
lect a unique periodic solution among those which 
satisfy ([7]), i.e. to fix a base point x^ = x^(0) on 
the cycle O^ (see Figure [2]). For this we require 
the end-point of the connection to belong to a 
plane orthogonal to the vector f{x~^{0),a), and 
the starting point of the connection to belong to 
a plane orthogonal to the vector /(a;~(0), a). This 
allows the base points x^(0) to move freely and 
independently upon each other along the corre- 
sponding cycles O^. 

3.3 The complete BVP 

The complete truncated BVP to be solved numer- 
ically consists of 

x^-T^f{x^,a) = 0, (12a 

x^{0)-x^{l) = 0, (12b 

w^ + T^f^{x^,a)w^ + \^w^ = 0, (12c 



w^il)-s^w^iO)-- 


= 0, 


(12d 


{w^{0),wH0))-l- 


= 0, 


(12e 


u-Tf{u,a) -- 


= 0, 


(12f 


(/(x+(0),«),«(l)-x+(0)) = 


= 0, 


(12g 


(/(x-(0),a),n(0)-x-(0)) = 


= 0, 


(12h 


{w+{0),u{l)-x+{0))-- 


= 0, 


(12i 


{w-iO),uiO)-x~iO))-- 


= 0, 


(12j 


u{0)-x~{0)\^-s^-- 


= 0, 


(12k 



where the last equation places the starting point 
u(0) of the connection at a small fixed distance 
e > from the base point a;~(0). The time vari- 
able is scaled to the unit interval [0, 1] , so that 
both the cycle periods T^ and the connection 
time T become parameters. Hence, besides a 
component of a, there are five more parameters 
available for continuation: the connection time T, 
the cycle periods T^, and the multipliers A^. 



4 Starting strategies 

The BVP described in the previous section are 
only usable when good initial starting data are 
available. Usually, such data are not present. 
Here we demonstrate how initial data can be gen- 
erated through a series of successive continua- 
tions in AUTO, a method referred to as homotopy 
method, first introduced by Doedel, Friedman and 
Monteiro (1994) for point-to-point problems and 
extended to point-to-cycle problems in Part I of 
this paper. 

4.1 Saddle cycles 

The easiest way to obtain the limit saddle cy- 
cles O^, first calculate a stable equilibrium using 
software like maple, matlab or mathematica. 
Then, using AUTO, continue this equilibrium up 
to an Andronov-Hopf bifurcation, where a stable 
limit cycle is generated. A continuation of this 
cycle can result in the detection of a fold bifurca- 
tion for the hmit cycle. This will yield a saddle 
limit cycle. 

4.2 Eigenfunctions 

In order to obtain an initial starting point for the 
connecting orbit we require knowledge about the 
unstable manifold of the saddle limit cycle 0~ . 
Also, we need the linearized adjoint "manifolds" 
to understand how the connecting orbit leaves 0~ 
and approaches O^ (or the same cycle in the ho- 
moclinic case). For this, we look at the eigendata. 
First consider the periodic BVP for 0~, 



X —T f\x , a) = 
x-(0)-x-(l) = 



(13) 



to which we add the standard integral phase con- 
dition 



{^oldi'^)^^ M) =0' 



(14) 



as well as a BVP similar to (IHl), namely 

V - T~ f^{x~ , a)v = 0, 

v{l) - fiv{0) = , 

{v{0),v{0))-h = 0. 



(15) 



In f lT^ . x~i^ is a reference periodic solution, e.g. 
from the preceding continuation step. The pa- 
rameter h in (IT^ is a homotopy parameter, that 
is set to zero initially. Then, (fT5l) has a trivial 
solution 

t;(t) = 0, h = 0, 

for any real /i. This family of the trivial solutions 
parametrized by n can be continued in auto us- 
ing a BVP consisting of (1131), dill), and (HSl) with 
free parameters (/i, h) and fixed a. The unsta- 
ble Floquet multiplier of 0~ then corresponds to 
a branch point at /i = /i~ along this trivial so- 
lution family. AUTO can accurately locate such 
a point and switch to the nontrivial branch that 
emanates from it. This secondary family is con- 
tinued in (yu, h) until the value h = 1 is reached, 
which gives a normalized eigenfunction v~ cor- 
responding to the multiplier /i~. Note that in 
this continuation the value of // remains constant, 
/i = /i~, up to numerical accuracy. For the initial 
starting point of the connection we use t'~(0). 

The same method is applicable to obtain the 
nontrivial scaled adjoint eigenfunctions w^ of the 
saddle cycles. For this, the BVP 



W^ _|_ J^ijTj^^^^ Q,j^ 



■,± 



r.± 



±.„± 



3±W±(0) 



w 



±/ 











(16) 

where s* = sign(/i^), replaces (IT^ . A branch 
point at X^ then corresponds to the adjoint mul- 
tiplier cr^e'^i . After branch switching the desired 
eigendata can be obtained. 

4.3 The connection 

Time-integration of ([1]), in matlab for instance, 
can yield an initial connecting orbit, however, this 
only applies for non-stiff systems. Nevertheless, 



mostly when starting sufficiently close to the ex- 
act connecting orbit in parameter space the method 
of successive continuation (Doedel, Friedman and 
Monteiro, 1994) can be used to obtain an initial 
connecting orbit. 

Let us introduce a BVP that is a modified 
version of (fT2|) 



X 



x^(l] 



w- 



~T^f{x^,a 
x±(0) • 

+ T^f^{x^,a)w^ + X^w^ 



w^{l) - s^w^{0) 



it — Tf{u, a) 

(/(x+(0),a),n(l)-x+(0))-/in 

(/(x-(0),a),M(0)-x-(0))-/ii2 

{w+iO),uil)-x+iO))-hi 

{w~{0),u{0)-x~{0))-h22 



(17a) 
(17b) 
(17c) 
(17d) 
(17e) 
(17f) 

(17g) 

(17h) 

(17i) 

(17j) 
(17k) 



where each <l>^ in 017cj) defines any phase condi- 
tion fixing the base point x^(0) on the cycle O^. 
An example of such a phase condition is 



$" 



X 



x,(0) 



A, J , 



where Qj is the jth-coordinate of the base point of 
0~^ at some given parameter values. Furthermore, 
hjk, j, k = 1,2, in (117hp - (117k[) are homotopy pa- 
rameters. 

For the approximate connecting orbit a small 
step e is made in the direction of the unstable 
eigenfunction v~ of the cycle 0~: 



{r)=x-{0) + ev-{0)e''^^'\ re [0,1], (18) 



u 



which provides an approximation to a solution of 
u = T~ f{u, a) in the unstable manifold W^ near 
O". After collection of the cycle-related data, 
eigendata and the time-integrated approximated 
orbit, a solution to the above BVP can be con- 
tinued in (T, /ill) and (T, /112) for fixed value of 
a in order to make hu = hi2 = 0, while u{l) 
is near the cycle O^, so that T becomes suffi- 
ciently large. In the next step, we then try to 



make /z2i = /i22 = 0, after which a good approxi- 
mate initial connecting orbit is obtained. 

This solution is now used to activate one of 
the system parameters, say ai, and to continue a 
solution to the primary BVP flT^ . Then, if neces- 
sary after having improved the connection first by 
a continuation in T, continuation in {ai,T) can 
be done to detect limit points, using the standard 
fold-detection facilities of auto. Subsequently a 
fold curve can be continued in two parameters, 
say (ai,Q;2), for fixed T using the standard fold- 
following facilities in auto. 

5 Implementation in AUTO 

Our algorithms have been implemented in auto, 
which solves the boundary value problems using 
super convergent orthogonal collocation with adap- 
tive meshes, auto can compute paths of solu- 
tions to boundary value problems with integral 
constraints and non-separated boundary condi- 
tions: 



f/(r)-F(f/(r),/3) 
b{U{0),U{l),(3) 



0, re [0,1], 
0, 



q{U{r),/3)dT = 



(19a) 
(19b) 

(19c) 



where 



[/(■), F{-, ■) G M"^ b{; ■) G W'^ g(-, ■) G W^% 
and 

p G M"■f^ 

as Ufp free parameters /? are allowed to vary, where 

nfp = Ubc + nic-nd + l . (20) 

The function q can also depend on F, the deriva- 
tive of U with respect to pseudo-arclength, and 
on U, the value of U at the previously computed 
point on the solution family. 

For our primary BVP problem (fT2!) in three 
dimensions we have 



Ud = 15, Ui 



0, 



and Ubc = 19, so that any 5 free parameters are 
allowed to vary. 



6 Example: food chain model 

In this section we describe the performance of the 
BVP-method for the detection and continuation 
of a cycle-to- cycle connecting orbit in the stan- 
dard food chain model, also used in Part I of this 
paper. 

6.1 The model 

The three-level food chain model from theoreti- 
cal biology, based on the Rosenzweig-MacArthur 
(1963) prey-predator model, is given by the fol- 
lowing equations 



Xi(l -Xi) - /l(xi,X2) , 
fl{Xl,X2) -diX2 - f2{x2,X3) 
/2 (3:2,3:3) -d2X3 , 




:2r 



with Holling Type-II functional responses 



and 



flixi,X2) 



f2{x2,X3) 



aiXiX2 
1 + biXi 



02X2X3 



1 + ^22:2 

This standard model has been studied by several 
authors, see e.g. Kuznetsov and Rinaldi (1996) 
and Kuznetsov et al. (2001) and references there. 

The death rates di and ^2 are often used as 
bifurcation parameters ai and 02, respectively, 
with the other parameters set at oi = 5, 02 = 0.1, 
hi = 3, and 62 = 2. For these parameter values 
the model displays chaotic behaviour in a given 
parameter range of di and (^2 (Hastings and Pow- 
ell, 1991; Klebanoff and Hastings, 1994; McCann 
and Yodzis, 1995). The region of chaos can be 
found starting from a fold bifurcation at for in- 
stance di ^ 0.2080452, ^2 = 0.0125, where two 
limit cycles appear. The stable branch then un- 
dergoes a cascade of period-doubhngs (see Fig- 
ure [3]) until a region of chaos is reached. 

Previous work by Boer et al. (1999, 2001) has 
shown that the parameter region where chaos oc- 
curs is intersected by homoclinic and heteroclinic 
global connections, and that this region is partly 





1 1 


,T;,„™(2x) 
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1 



0.2 



0.22 



0.24 



0.26 



0.28 



0.3 



di 



Figure 3: One-parameter bifurcation diagram for 
d2 = 0.0125. The equilibrium is indicated as x^. 
The dashed hne Lp is the Xs-value of the local 
minimum of an unstable (saddle) limit cycle. At 
the point T^ this limit cycle coincides with a sta- 
ble limit cycle, of which the local minimum of X3 
is also shown. The stable limit cycle undergoes 
period doublings until chaos (the dense regions) 
is reached. The two dense regions are separated 
by a region where homoclinic cycle-to-cycle con- 
nections to the limit cycle Lp exist. Both chaotic 
regions are bounded by a limit point of the ho- 
moclinic connection, indicated by Thom- Observe 
that near the right chaotic region, between two 
limit points exist secondary connecting orbits to 
the cycle. One of these limit points coincides with 
the limit point of the primary connecting orbit 
that forms the boundary of the region of chaos 
(hence Thom twice); after Boer et al., 1999). 



x. 



12.5 



11.5 



10.5 




Figure 4: Phase plot of the homoclinic cycle-to- 
cycle connection. 



bounded by a homoclinic cycle-to-cycle connec- 
tion, as shown in Figure [31 These results were 
obtained numerically using multiple shooting. 

6.2 Homotopy 

Using the technique discussed in this paper we 
first find the saddle limit cycle for di = 0.25, ^2 = 
0.0125. Since the cycle O is both O"*" and O", we 
use the same initial base point 

x^{0) = (0.839783,0.125284,10.55288) 

and the period T^ = 24.28225. The logarithms 
of the nontrivial adjoint multipliers are 



A^ 



-0.4399607 , A^ = 6.414681 



The starting point of the initial "connecting" or- 
bit is calculated by taking the base point x~(0) 
and multiplying the eigenfunction f~(0) hj e = 
-0.001 

u{0) = X- (0) + eir (0) , (22) 

where 

v~{0) = (-1.5855 ■10^^ 2.6935 -10"^ -0.99951) , 

and the resulting 

m(0) = (0.839789,0.125274,10.55324) . 

The connection time T = 503.168. 

To obtain a good initial connection we con- 
sider a BVP like fITTj) . with 6 free parameters: /i^, 
T^, T, and, in turn, one of the four homotopy pa- 
rameters /ill, hi2, h2i, /i22- The selected boundary 
conditions 017cp are 



$+[x] =a;^(0)- 0.125274 



and 



$- 



X 



x^(0)- 0.839789 



so, the first condition uses the a;2-coordinate of 
the initial base point selected on the cycle, while 
the second condition uses the a;i-coordinate of the 
initial base point. Observe, that this selection is 



11.5 - 




10.5 




Figure 5: Profiles of tlie liomoclinic cycle-to-cycle connections for d2 = 0.0125. The left panel 
compares the two profiles of the connections for di = 0.2809078, the right panel compares the four 
profiles for di = 0.27850, which is between the limit points for the primary and secondary branches. 
The connection time T is scaled to one. 



somewhat arbitrary and that one can select other 
base point coordinates. 

In the continuation we want hu = hu = 
^21 = ^22 = 0. However, there are several solu- 
tions, that correspond to connecting orbits with 
different numbers of excursions near the limit cy- 
cle, both at the starting and the end-part of the 
orbit. Observe that the success of the future con- 
tinuation in (^1,(^2) seems to depend highly on 
the number of excursions near the cycle at the 
end-point of the connecting orbit. In the food 
chain model a decrease in ^2 is accompanied by 
a decrease in the numbers of excursions near the 
cycle at the end-point of the connection, like a 
wire around a reel. If this number is too low, a 
one-parameter continuation in di 2 will yield in- 
correct limit points. Also, two-parameter contin- 
uations in ((^1,^2) will most likely terminate at 
some point. Hence a starting orbit is selected 
with a sufficient number of excursions near the 
cycle at the end-point, with T = 454.04705 and 
e^ = 0.069414. 

6.3 Continuation 

The continuation of the connecting orbit can be 
done in di 2 using the primary BVP flT21) . Equa- 
tion (112kp ensures that the base points x^(0) G 



O"^ become different (and so do the periodic so- 
lutions a:^(t)). 

First, however, using this BVP, the connec- 
tion can be improved by increasing the connec- 
tion time, for the same reason as mentioned above 
with regard to the number of excursions near the 
cycle at the end-point. The increase in T re- 
sults in an increase of the number of excursions 
near the cycle at the end-point of the connecting 
orbit. Then, the continuation in di^2,T results 
in the detection of four limit points, of which 
two are identical. Figure [3] shows a continua- 
tion in di for fixed ^2 = 0.0125, that detects 
four limit points in one "run". Observe that in 
this way not only the primary [di = 0.2809078, 
twice, and di = 0.2305987), but also the sec- 
ondary [di = 0.2776909) branch is detected. Fig- 
ure E] shows the profiles of the connecting orbits 
for d2 = 0.0125. Observe that for the region of 
0.2776909 < di < 0.2809078 there are four differ- 
ent connecting orbits with the same connection 
time T (see right panel). 

Using the standard fold-following facilities for 
BVPs in AUTO, both critical homoclinic orbits 
can be continued in two parameters [di, ^2)- Along 
these orbits the stable and unstable invariant man- 
ifold of the cycle are tangent. Starting from di = 
0.2809078 we continue the primary branch. The 
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Figure 6: Two-parameter curves of the primary (left) and secondary (right) homochnic cycle-to- 
cycle connections of the food chain model. Depicted is also the fold bifurcation curve of a limit cycle 
(dashed). 



secondary branch is continued from d\ 
Both curves are depicted in Figure [61 



0.2776909. 



7 Discussion 

Our continuation method for cycle-to-cycle con- 
nections, using homotopies in a boundary value 
setting, is a modified method proposed in our 
previous paper for the continuation of point-to- 
cycle connections (Doedel et al., 2007). The re- 
sults discussed here seem to be both robust and 
time-efficient. Detailed auto demos performing 
the computations described in Section 6 are freely 
downloadable from 

www . bio . vu . nl/thb/research/pro j ect/globif 
Provided that the cycle has one simple unsta- 
ble multiplier, the proposed method can be ex- 
tended directly to homoclinic cycle-to-cycle con- 
nections in ra-dimensional systems. 
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